R Notebook

library(psych)
library(corrplot)
library(dplyr)
library(tidyverse)
library(stringr)
library(base)
library(sf)
library(scales)
library(vars)
library(lattice)
library(urca)
library(ggplot2)
library(MASS)
library(haven)
library(joineR)

Loading required package: survival

library(pROC)

Type 'citation("pROC")' for a citation.

Attaching package: ‘pROC’

The following objects are masked from ‘package:stats’:

    cov, smooth, var

library(MESS)

Attaching package: ‘MESS’

The following object is masked from ‘package:pROC’:

    auc

library(arules)

Loading required package: Matrix

Attaching package: ‘Matrix’

The following objects are masked from ‘package:tidyr’:

    expand, pack, unpack


Attaching package: ‘arules’

The following object is masked from ‘package:dplyr’:

    recode

The following objects are masked from ‘package:base’:

    abbreviate, write

read file

data=read.csv('new_data.csv')

Data exploration

There are 12 variables in the data with 485 observations. Including 1 date variable, 2 response variables, and 9 predictors. All variables are numeric type except variale ‘date’. 1. date 2. new_cases: the daily new confirmed COVID cases in Netherlands 3. new_deaths: the daily new deaths of COVID patients in Netherlands 4. ICU_patients: the COVID patients admitted to ICU per day in hospital 5. Hosp_patients: the COVID patients admitted to hosptial per day 6. stringency_index: 7. temp: daily temperature 8. humidity: daily humidity 9. precip: the amount of rains per day 10. windspeed: daily windspeed 11. transit_stations: 12: residential:

问题： transit 都是负数，需不需要调整从0开始 判断人们乘坐交通的变化

checking missingness

sum(is.na(data$hosp_patients))

[1] 9

There are 18 missing values in the data, locate in ICU patients and Hospital patients columns. (The government website does not have COVID records for some days during 2021- 2022)

correlations

For hospital admission 1.hospital admission and humidity, hospital admission and precipitation have small correlation less than 0.1 2.hospital admission and new deaths, hospital admission and resident are highly positive correlated 3.hospital admission and temperature, hospital admission and transit are highly negative correlated

for hospital admission: fit var model including all variables, comparing to the var model without humidity, precipitation, and wind speed.

For ICU admission 1.ICU admission and humidity, ICU admission and precipitation, ICU and wind speed have small correlation less than 0.1. 2.ICU admission and new deaths, ICU admission and stringency, ICU admission and resident are highly positive correlated 3.ICU admission and temperature, ICU admission and transit are highly negative correlated

for ICU admission: fit var model including all variables, comparing to the var model without humidity, precipitation and wind speed.

new deaths highly positive correlated to stringency_index, and residential. Negative correlated to temperature and public transition.

strigency_index increases, residential increases, transition decreases. (temperature has opposite relation)

transit and residential as highly negative correlated.

Plots

plot(data$new_cases) # 2022-02-07 has 380399 cases, (outlier, remove and graph again) volatity

plot(data$new_cases[-403])

plot(data$new_deaths)

# more deaths at the beginning 
plot(data$icu_patients)

plot(data$hosp_patients)

plot(data$stringency_index)

plot(data$temp)

plot(data$humidity)

plot(data$precip)

plot(data$windspeed)

plot(data$transit_stations)

plot(data$residential)

# ICU, hosp, stringency, residential trends are similar, transition has opposite trend.

stringency 会导致lockdown， residential的变化和stringency相似， 从而影响hosp admission。

summer low COVID, winter high COVID 湿度看着和hosp admission， icu走势有点像， 风速和降雨量看不出来。

distribution

hist(data$icu_patients)
hist(data$hosp_patients)
hist(data$new_cases)
hist(data$new_deaths)
hist(data$stringency_index)
hist(data$temp)
hist(data$humidity)
hist(data$precip)
hist(data$windspeed)
hist(data$transit_stations)
hist(data$residential)

heteroscedasticity

lm1=lm(hosp_patients~., data=data[,c(2,3,5,6,7,8,9,10,11,12)])
summary(lm1)
# number of new deaths, temperature, humidity, and residential have significant impact on hospital admission.
# new deaths increases, more hospital admissions
# temperature and humidity get higher, the less hospital admission
# hospital admission increases as residential increases.

plot(lm1)
# heteroscedasticity detected, maybe garch model can be used in time series model (further hetero test need to be done)

如果将温度， 湿度， 风速降雨量等调整为categorical数据， 是否能观察出相关性

data_cat=data
data_cat$temp_cat_int <- discretize(data_cat$temp, method="interval", breaks=4)
data_cat$temp_cat_fre <- discretize(data_cat$temp, method="frequency", breaks=4)
data_cat$temp_cat_clu <- discretize(data_cat$temp, method="cluster", breaks=4) #KNN
summary(data_cat$temp_cat_int)
summary(data_cat$temp_cat_fre)
summary(data_cat$temp_cat_clu)
# based on discretization, we used 3 different method to categorize the temperature into 5 level, select the result from KNN model
# [-6,4): cold
# [4,9)：cool
# [9,14)： warm
# [14,23]： hot

data_cat$wind_cat_int <- discretize(data_cat$windspeed, method="interval", breaks=4)
data_cat$wind_cat_fre <- discretize(data_cat$windspeed, method="frequency", breaks=4)
data_cat$wind_cat_clu <- discretize(data_cat$windspeed, method="cluster", breaks=4) #KNN
summary(data_cat$wind_cat_int)
summary(data_cat$wind_cat_fre)
summary(data_cat$wind_cat_clu)
#[7.6,19.5):low
#[19.5,29): mid
#[29,43.4): strong
#[43.4,82.5]: very strong

data_cat$pre_cat_int <- discretize(data_cat$precip, method="interval", breaks=4)
data_cat$pre_cat_fre <- discretize(data_cat$precip, method="frequency", breaks=4)
data_cat$pre_cat_clu <- discretize(data_cat$precip, method="cluster", breaks=4) #KNN
summary(data_cat$pre_cat_int)
summary(data_cat$pre_cat_fre)
summary(data_cat$pre_cat_clu)

data_cat$hum_cat_int <- discretize(data_cat$humidity, method="interval", breaks=4)
data_cat$hum_cat_fre <- discretize(data_cat$humidity, method="frequency", breaks=4)
data_cat$hum_cat_clu <- discretize(data_cat$humidity, method="cluster", breaks=4) #KNN
summary(data_cat$hum_cat_int)
summary(data_cat$hum_cat_fre)
summary(data_cat$hum_cat_clu)
#[46.6,68.4): level 1
#[68.4,79.8): level 2
#[79.8,88.3): level 3
#[88.3,98.8]: level 4

convert the numerical variables into categorical variables

data_cat$temp_cat=ifelse(data$temp<4, 'cold', 
                     ifelse(data$temp<9, 'cool',
                            ifelse(data$temp<14, 'warm', 'hot')))
data_cat$wind_cat=ifelse(data$windspeed<19.5, 'low', 
                     ifelse(data$windspeed<29, 'mid',
                            ifelse(data$windspeed<43.5, 'high', 'strong')))

plot the categorical variables vs number of hospital patients

fit a lm using categorical data

lm2=lm(hosp_patients~., data=data_cat[, c(2,3,5,6,11,12, 25, 26)])
summary(lm2)
# when temperature is higher than 9 degrees, it has an effect on the COVID hospital admission.
# wind does not have an impact on hospital admission